rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
#9.1 生存对象
library(survival)
data(ovarian)
summary(ovarian)
?ovarian
#因子转换
ovarian$resid.ds <- factor(ovarian$resid.ds,
                           levels = c(1,2),
                           labels = c("no", "yes"))
ovarian$rx <- factor(ovarian$rx,
                     levels = c(1, 2),
                     labels = c("A", "B"))
ovarian$ecog.ps <- factor(ovarian$ecog.ps,
                          levels = c("1", "2"),
                          labels = c("good", "bad"))
#查看连续型变量age的分布
hist(ovarian$age)
#图中age非对称；考虑到结果的易解释性，
#使用cut将age进行因子转换
ovarian$agegr <- cut(ovarian$age, 
                     breaks = c(0, 50, 75),
                     labels = c("<=50", ">50"))
table(ovarian$agegr)
#用函数Surv创建生存对象
surv.obj <- Surv(time = ovarian$futime, event = ovarian$fustat)
surv.obj

#9.2 生存率的估计与生存曲线
#乘积极限法
#生存率的Kaplan-Meier估计
#~1表示没有自变量
survfit(surv.obj~1)
#因为该数据集的样本量较小，生存曲线并没有穿过生存率在50%对应的水平线
#无法计算中位时间置信区间的上限
surv.all <- survfit(surv.obj~1)
summary(surv.all)
#删失数据并没有显示;显示删失数据
summary(surv.all, censored = TRUE)
#生存曲线绘制
#mark.time = TRUE，为删失数据添加标记“+”
#虚线为生存率的置信区间
#不想显示置信区间，conf.int=FALSE
plot(surv.all, mark.time = TRUE, col = "purple")


#9.3 生存率的比较
#通过在函数survit的公式中添加一个因子变量
#来获得该变量的不同水平下的生存信息
#例如，欲获得不同治疗方法下生存率的估计
surv.treat <- survfit(surv.obj~rx, data = ovarian)
summary(surv.treat)
#在同一图中显示多条生存曲线
plot(surv.treat, mark.time = TRUE, conf.int = TRUE,
     lty = c(1,2), col = c("blue","red"))
#install.packages("survminer")
library(survminer)
ggsurvplot(surv.treat, data = ovarian, pval = TRUE)

#以上结果表明，B法生存率高于A法，
#这种差异是偶然的还是由治疗方法的不同引起，需要进行统计学检验
#最常用的是时序检验（log rank test）
#基本思想：先计算出不同时间，两种治疗方法的暴露人数和死亡人数
#在两种治疗方法效果相同的假设下计算出期望死亡人数
#如果不拒绝零假设，则实际观测值和期望值的差异不会很大
#如果差异过大则不能认为是由随机误差引起的差异
#用卡方检验做判断；时序检验可以用函数survdiff来实现
survdiff(surv.obj~rx, data = ovarian)

#由p=0.3可知，两种治疗方法的生存率差异没有统计学意义
#但是组中被研究的对象存在其他可能影响生存结果的因素
#例如年龄、疾病残留情况、患者的ECOG评分
#此时需要用回归分析，把这些因素作为协变量来解释两组生产结果的差异
#在平衡（调整）了这些潜在的影响因素后，
#两组生存时间的分布的比较才能更加准确

#9.4 Cox回归
#所有的参数回归模型需要对风险函数做出假设，
#对每个时间点的生存概率进行估计
#Cox回归没有关于风险函数的假设
#有个函数、公式，自己去看

#9.4.1 建立Cox回归模型
#将所有的协变量都包含进来建立Cox回归模型
cox1 <- coxph(surv.obj~rx+resid.ds+agegr+ecog.ps,
              data = ovarian)
summary(cox1)
#变量选择
drop1(cox1, test = "Chisq")
#ecog.ps对应的p值最大，移除
cox2 <- coxph(surv.obj~rx+resid.ds+agegr,
              data = ovarian)

step.cox1 <- step(cox1)
#结果表明，模型cox2就是最优模型

#9.4.2 比例风险假定的检验

#验证Cox回归的比例风险假定
cox.zph(cox2)
#各个自变量的检验和全局检验的p值都大于0.05，没有违背比例风险的假定
summary(cox2)
#在调整了疾病残留情况和年龄的情况下，
#两种治疗方法下的死亡风险的差异由统计学意义
#不同疾病残留情况以及不同年龄组之间的死亡风险的差异都没有统计学意义

#9.4.3 生存的预测
#建立新数据
newdata <- data.frame(rx = c("A","B"),
                      resid.ds = c("no","no"),
                      agegr = c(">50",">50"))
newdata
#预测
hr <- predict(cox2, newdata = newdata, type = "risk")
hr
hr[1]/hr[2]
#表明年龄大于50岁的且没有疾病残留的患者，
#采用治疗方法A的死亡风险是B的3.599421倍

#生存曲线拟合
cox.fit <- survfit(cox2, newdata = newdata, type = "kaplan-meier")
plot(cox.fit, lty = c(1,2), col = c(2,4))
title(main = "Cox survival curve by treatment
      for age > 50, no residual disease patients",
      xlab = "Duration in days",
      ylab = "Survival probability",
      las = 1)
#10为离ylab的距离，0.3为高度
legend(5,0.3,c("Treatment A","Treatment B"),
       lty = c(1,2), col = c(2,4))

